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Abstract 

For accurate simulations of rarefied gas flows around moving obstacles, we propose a 
cut cell method on Cartesian grids; it allows exact conservation and accurate treatment 
of boundary conditions. Our approach is designed to treat Cartesian cells and various 
kind of cut cells by the same algorithm, with no need to identify the specific shape 
of each cut cell. This makes the implementation quite simple, and allows a direct 
extension to 3D problems. Such simulations are also made possible by using an adaptive 
mesh refinement technique and a hybrid parallel implementation. This is illustrated 
by several test cases, including a 3D unsteady simulation of the Crookes radiometer. 

Keywords: kinetic equations, deterministic method, immersed boundaries, cut cell method, 
rarefied gas dynamics 


1 Introduction 

In gas dynamic problems, the rarefied regime appears when the mean free path of the 
molecules of the gas is of the same order of magnitude as a characteristic macroscopic 
length. The flow has to be modeled by the Boltzmann equation of the kinetic theory of 
gases. Most of numerical simulations for rarefied flows are made with the stochastic DSMC 
method [6], especially for aerodynamical flows in re-entry problems. In the past few years, 
several deterministic solvers have been proposed, that are based on discretizations of the 
Boltzmann equation or simplified models, like BGK, ES-BGK, or Shakhov models [2H]- 
They are efficient for accurate simulations, multi-scale problems, or transitional flows, for 
instance. 

A recent issue is the account of solid boundary motion in rarefied flow simulations. 
This is necessary to simulate flows around moving parts of micro-electromechanical systems 
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(MEMS) [IZIES], as well as flows inside vacuum pumps. A fascinating illustration of rarefied 
flows with moving boundaries is the Crookes radiometer, subject of many debates from the 
late 19th to early 20th century [21]. Recent deterministic simulations help to understand 
the origin of the radiometric forces [3S111I1I121IS]. The numerical simulation of the Crookes 
radiometer is difficult because the motion of the vanes is induced by gas/solid interaction 
(like thermal creep), which means that an accurate prediction of the flow in the vicinity of 
the boundary is needed in order to predict the correct velocity of the vanes. 

There are several numerical methods for moving boundary problems designed for com¬ 
putational fluid dynamics: some of them have recently been extended to deterministic dis¬ 
cretizations of kinetic models, and can be divided in two main categories. 

First, with body htted methods, the mesh is adapted at each time step so that the bound¬ 
ary of the computational domain always £t with the physical boundary: moving mesh |13] 
and ALE methods [IHl HH] fah into this category. Despite their extensive use in computa¬ 
tional fluid dynamics, very few similar works have been reported in kinetic theory, except by 
Chen et ah [in|. Methods of the second category are based on Cartesian grid computations 
and are usually referred to as immersed boundary methods m- The mesh does not change 
during computations, and hence does not £t with the physical boundary. Special treatment 
is applied on mesh cells that are located close to the boundary in order to take its motion 
into account. Various extensions of these methods to kinetic theory have been proposed by 
several authors in [21 |3ll HU 0]. Two recent variants are the inverse Lax-Wendroff immersed 
boundary method proposed by Filbet and Yang [I6] and the Cartesian grid-based unihed gas 
kinetic scheme of Chen and Xu [8|: the boundary motion is not taken into account in these 
two works, but these methods could in principle be extended to this kind of problem. We 
also mention the Lagrangian method: while it falls into the first category in CFD, it does 
not in kinetic theory. Indeed, whatever the motion of the mesh, the distribution function 
has to be interpolated at the foot of the characteristic for each microscopic velocity. The 
accuracy of these methods have been shown in [35lll^ for one dimensional problems. Finally, 
we mention that moving boundary flows can also be treated with DSMC solvers: see, for 
instance, [201 ESI EOl SO] ■ 

In this paper, we try to mix the advantages of body htted and Cartesian methods: we 
present a cut cell method for computing rarehed gas hows around moving obstacles. The 
cut cell method belongs to the Cartesian grid based methods and has been widely used in 
computational huid dynamics [22]. However, this is the hrst extension to moving boundary 
problems in kinetic theory (complex 3D stationary DSMC simulations have already been 
investigated in [221 SB])- This approach is well suited to deterministic approximations of 
the Boltzmann equation and is easy to implement because of the Cartesian structure of 
the mesh. Moreover, this is, up to our knowledge, the only immersed boundary method to 
be conservative. The versatility and robustness of the technique is illustrated by various 
2D hows, and by the simulation of the unsteady rotation of the vanes of a 3D Crookes 
radiometer. This article is an extended version of our work announced in [IT]. Here, the 
Boltzmann collision operator is replaced by BGK like models, which is approximated by a 
discrete velocity method. However, this is not a restriction: other collision operators could 
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be used, and any velocity approximation (like the spectral method) could be used. 

Generally, the problem of cut cell methods is that it is difficult to take into account the 
various shapes of cells that are cut by the solid boundary: for instance, in 2D, a cut cell 
can be a triangle, a quadrangle, or a pentagon, and this is worse in 3D. Here, we propose a 
simple representation of these cells by using the notion of virtual cells that are polygons (or 
polyedrals) with possibly degenerated edges (or faces). This makes the treatment of any cut 
cell completely generic: in the implementation, the different kinds of cut cells and the non 
cut cells are treated by the same algorithm. This makes the extension of the method to 3D 
problems very easy. However, to make large scale 3D simulations possible, we also use an 
adaptive mesh rehnement (AMR) technique and a special parallel implementation. 

The outline of our paper is as follows. In section we give the governing equations of 
rarehed gas flows and introduce some notations. Our cut cell method is presented in section]^ 
for 2D problems. It is validated on three different numerical examples in section Then, in 
section our algorithm is extended to 3D simulations, and a 3D unsteady simulation of the 
Crookes radiometer is presented. Finally, some conclusions and perspectives are discussed 
in section Technical details like computations of geometric parameters of the cells are 
presented in the Appendix. 


2 Rarefied gas dynamics 


2.1 Boltzmann equation 

In rarehed regimes, a monoatomic gas is described by the Boltzmann equation: 

^ + v-VF = Q(F). ( 1 ) 


The distribution function F(f, T, v) is the mass density of molecules at time t that are located 
at the space coordinate T G and that have a velocity u G For our approach, it is more 
relevant to look at the integral form of ([^ in a time dependent volume V{t). The Reynolds 
transport theorem leads to: 


m 


FdV + 


'V 



■nFdS 


[ Q{F) dy 
Jv 


( 2 ) 


where dV{t) is the surface of the volume V{t). Let T be a point of this surface: it is moving 
at a velocity w(t,x) and the vector n(t,x) is the outward normal vector to the surface at 
this point. 

The density p, momentum pu, total energy E, stress tensor S and heat hux q, are 
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computed by the first moments of the distribution function with respect to the velocity: 



F(t, X, v) dvxdvydvz, 


S 



(n — u) ® {v 


u) F(t, X, v) dvxdvydvz, 


q = 


-{v - u)\\v - u\\‘^ F{t,x,v) dVxdVydVz, 


( 3 ) 


where the norm is dehned by ||n|p = F v'^ + v^. The temperature T of the gas is related 
to the the energy by the relation E = IpHnp + ^pRT, where R is the gas constant dehned 
as the ratio between the Boltzmann constant and the molecular mass of the gas. Moreover, 
the pressure is computed with the standard equation of state for ideal gases: P = pRT. 

When a gas is at rest, which means in equilibrium state, the molecules are uniformly 
distributed around the macroscopic velocity and the distribution function is a Gaussian 
function called Maxwellian: 

Mip,u,nv) = 

The Boltzmann equation implies that the total variation of F comes from the collisions 
between particles, that are modeled by the Boltzmann operator Q{F). This operator is 
computationally expensive and several models have been introduced to make its computa¬ 
tion easier. Simplest models consists in a relaxation of the distribution function towards a 
corresponding equilibrium function S: 

Q(F) = 1(£ - F), 

T 

where r is a relaxation time (see various dehnitions in section]^. Bhatnagar, Gross and 
Krook |5] proposed to take £ equal to the local equilibrium state, that is 


£ = M[p,u,T]. 


The corresponding BGK model conserves mass, momentum and total energy since the hrst 
three moments of the Maxwellian are the same as those of the distribution function. A 
Ghapman-Enskog expansion [7] gives relations between viscosity p, heat conduction k and 
relaxation time. In this case, the expansion yields p = r P and k = ^Rr P. It may be seen 
that the BGK model necessarily leads to Prandtl number Pr = ^pR /k equal to 1. However 
for most of gases, it is physically found to be less than this. More complex functions £ 
such that in Holway [20] and Shakhov [2H| models have been developed to obtain the correct 
Prandtl number. For instance, for the Shakhov model, £ reads 


£ = M[p,u,T] 


1 -I- (1 — Pr)(u — u) ■ q 


^-5)/(5Pisr) 
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In this article, solid wall interactions are taken into acconnt by the standard fnlly diffuse 
reflection. This model states that all particles that collide with a boundary are absorbed by 
the wall and re emitted with a Maxwellian distribution: 


e r,n e Vin) = (5) 

where and are the temperature and the velocity of the boundary T at position x. The 
coefficient 0 is computed in order to set the net mass flux across the wall to zero: 




n. 


' l^GVout 


^ Fdv^dvydv^ 


Iv&Vi, 


{v - u^) ■n^M[l, Uw,Tj\ dv^dvydv, 


The set of incoming velocities is defined by Vm = {w such that {y — u^) ■ Hw < 0}, where 
fiy^ is the normal vector to the boundary, pointed to the wall. Similarly, the set of outgoing 
velocities is Vout = {h* such that (n — Uu}) ■ > 0}. Note that the boundary condition is 

defined only for the relative incoming microscopic velocities. 


2.2 Reduced model 

For plane flows, the computational complexity of the Boltzmann equation can be decreased 
by the use of a standard reduced distribution technique HU This classical method has been 
extensively used for numerical computations of BGK and Shakhov models. First note that 
in plane flows, the third component of the macroscopic velocity is equal to zero, as well 
as Qz-, ^xz, and T^yz- From now on, we define the two dimensional variables 

^={x,y), v = {vx,Vy), u={ux,Uy), q={qx,qy), F = ( 

\^yx ^yy 

Let / and g be the reduced distribution functions defined by 

/(v) = [ F{v)dvz, ^(v) = [ ^vlF{v)dVz. 

Jr Jr ^ 

The macroscopic quantities can be computed from / and g. Indeed, the set of equations (|^ 
readily becomes 



p 

r 

1 

r 

■ 0 ■ 


pu 

= / 

V 

f{-v)dvxdvy+ / 

0 

g{-v) dVxdVy, 

E 

JR2 


Jr2 

1 


S 

q 



u) ® (v - u)/(v) dvxdvy, 

u) Q||v - u|p/(v) + ^(v)^ dVxdVy, 


( 6 ) 
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where ||v|p = + Vy. Multiplying by (1, \v‘l) the Boltzmann equation (|^ and then inte¬ 

grating the result with respect to Vz gives the following set of equations: 


4 //d^+ / (v-w)-n/d/= [-{F-f)dS, 

Js Jds Js ^ 


(7) 


In this case, dS{t) is the contour of the surface S{t). Each point x G dS{t) is moving at a 
velocity w{t,x) and n{t,x) is the outward normal vector to the contour at this point. The 
reduced equilibrium functions are dehned by 

RT 

f = M[p,u,T] and G = u, T], (8) 

for the BGK model and by 


F = M[p,u,T] 


1 4- (1 — Pr)(v — u) ■ q 


V — u 


RR 

G = —M[p,u,T] 


1 -|- (1 — Pr)(v — u) ■ q 


RT 

V — u 


-4j /{5PRT) 
2 


RT 


2 / {5PRT) 


(9) 


for the Shakhov model, where M[p, u, T] is the reduced Maxwellian given by 

12 ' 


M[p,u,T](v) = 


P 


27iRT 


exp 


V — u 


2RT 


To close this section, the boundary conditions (§ are written with the reduced distribution 
functions as: 

/(t,x G P, V G Vin) = (^M[1,U^,T^], 

RT (10) 

g{t,yie P, V G Vin) = 0—M[l,Un„Tn,], 


where 0 is computed by 


u„ 


n. 


, / dvz;dv^ 


' veVout 


u„,) • M[1, TJ\ dv^dv^^ 


'veVn 


In this formula, and Tyj are the velocity and temperature of the point x that belongs 
to the boundary P, and n^, is the normal vector to the boundary pointed to the wall. 
The relative incoming and outgoing velocities at this point are therefore dehned by Vin = 
{v such that (v — u^,) ■ < 0} and Vout = {v such that (v — u^,) • n^, > 0}. 


3 The cut-cell method for two dimensional problems 


In this section, we present a numerical method to simulate plane hows with moving bound- 

The discretization of each variable 


aries. 


The governing equations are detailed in section 2.2 


(velocity, space, and time) is presented in separate sections. 
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3.1 Discrete velocity approximation 

The velocity space is discretized by a Cartesian grid. Let Vmin G and v^ax ^ be 
the lower-left and upper-right corners of this grid. The number of discrete velocities is 
N^, the velocity step is denoted by {Avx,Avy) = (v^ax — ^ram)/Af, and the velocity is 
Vp = Vmin + iPiAvx,P 2 Avy), such that p = P 2 N + Pi for all (pi,P 2 ) e [0, iV - 1]^. The 
approximation of the distribution function is dehned by fp{t,x.) = The set of 

equations Q is discretized with respect to v by the following set of 2N‘^ equations: 

4 / /p + / (vp-w)-nfpdl = [-{Fp- fp) d^, 

Js JdS Js 

— / 5'pdS'-F / (vp - w) ■ n^ipd/= / -(Cp - ^fp) dS*. 

Js Jas Js T 

The macroscopic quantities are computed with ([^, where the integrals over are approx¬ 
imated by a sum over the discrete velocity points. They are therefore given by 


p 

N^-l 

1 

N^-1 

■ 0 ■ 

pu 

= E 

Vp 

fp ^Vx^Vy ^ ^ 

0 

E 

p=0 

-1 Iv 1 P 

L J 

p=0 

1 


9p AvxAvy, 


S 

q 


N^-l 

p=0 

N^-l 

p=0 


u) ® (vp - u)/p AvxAvy, 

u) f ^l|vp - u||Vp + fi'p j AvxAvy. 


( 12 ) 


Finally, the equilibrium functions Fp and Gp are computed either with ([^ or with 
However, the triplet (p, u, T) used in these formulas is obtained with a Newton algorithm 
that preserves the discrete conservation of Boltzmann equation, rather than with direct 
computation (12) of the macroscopic quantities. Note that instead of using the algorithm 
of [2Z] which is based on entropic variables, we use the algorithm of Titarev 


3.2 Space discretization 

3.2.1 Cartesian grid and cut cells 

Let n = [xmin, a^max] X [pmin, 2/max] denote the space computational domain. It is discretized 
by a Cartesian grid of {Nx -f- 1) x {Ny + 1) points. Their coordinates are computed for 
all {i,j) e [0,Nx] X [0,iVp] by = Xmm + (iAxjAy), where {Ax, Ay) = ([xmax - 

X„iin]/Nx, [i/max “ Uminl/Ny) aud x^in = (a^min, 2/min) • The Computational mesh is therefore 
made up by Nx x Ny rectangular cells: each cell is denoted by flij and its center is x,j. 

Since the computational domain is rectangular, physical boundaries do not necessarily 
£t with the mesh boundary. In order to simulate arbitrary shaped objects, solid and gaseous 
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domains fls(t) and fig(t) are introduced. They correspond to the solid and gaseous parts of 
the computational domain and hence = Qs{t) UQg{t). We point out that while and Qg 
are time dependent, is not. At time t > 0, a rectangular cell can be in one of these 
three different states only: 

• Qij is a gas cell if it is completely contained in the gaseous domain: Qij fl Qs{t) = 0- 

• Qij is a solid cell if it is completely contained in the solid domain: Qij fl ^g{t) = 0. 

• Qij is a cut cell if it is partially contained in the gaseous domain and partially contained 
in the solid domain: flij fl fls(t) ^ 0 and flij fl flg(t) ^ 0. 

These three states of cells are shown in hgure 


3.2.2 Virtual cells 

To each cell is now associated a virtual cell Qij{t), which is the section of Qij contained 
in the gaseous domain: this reads fli,j(t) = flij nflgit). Whatever the state of the cell (that 
is to say gas, solid or cut), it is dehned with hve virtual edges, whose lenghts can be zero. 
Four of them, that are denoted by L-^i j(t) and £t with the lines of the Cartesian 

mesh. The last one, denoted by is a linear approximation of the solid boundary. If 

the virtual cell has less than five real edges, then at least one length is zero. Finally, we 
denote by Sij the area of the virtual cell, and \L\ will denote the length of any edge L. 

At a given time t” = nAt, all these parameters are denoted as follows: 




L” 1 . 


= L 


i±n,j 


A 


L ^.,1 

*J±2 


= L 


i,j± 


in Lig = mn n = s^,{e)■ 


Note that a difficult problem in the cut cell method is that cut cells can take many 
different shapes (mainly in 3D), which can make the code very complex. A key element of 
our approach is that all these different shapes are treated generically by using this notion 
of virtual cell with its 5 virtual edges. Indeed, all the cells are treated in the same way, 
whatever their state (gas, solid, cut cell) or shape. The different parameters of the three cell 
states are summarized below, and we refer to figure for three examples of cut cells: 


• gas cell: 

• solid cell: VL, = 0, i -I 




= 0 , 


= \n\ = ^n,g = o, 


AxAy, 


• cut cell: is a part of Qij, and the five lengths of the virtual cell ^.|, |LT^i|, 

and \L^j\ can take any value between 0 and Ay, Ax, and AAx'^ + Ay'^, respectively 
(see figure . 

All these parameters are computed with a levelset technique, as explained in appendix [A| 



3.2.3 Control volumes 


The notion of control volume is essential to avoid the use of very small virtual cells that 
would lead to prohibitively small time steps. The idea is to merge small virtual cells with 
larger neighboring cells when their areas are smaller than half of the area of a Cartesian cell. 

The control volume is constructed by recursion: we look at a given cut cell VLij whose 
center is inside the solid domain. The corresponding virtual cell j necessarily has an area 
smaller than ^AxAy, and it has to be merged with one of its non solid neighboring cells. 
This cell is chosen by looking at the largest non solid edge of 0^^: the neighboring cell 
that shares the same edge is chosen for merging (for instance in figure 3, top). If the 

corresponding neighboring cell has its center inside the gas domain, then the algorithm is 
stopped and the resulting control volume contains two virtual cells. It happens sometimes 
that the neighboring virtual cell is also too small (its center is inside the solid domain too): in 
this case, the same algorithm is used recursively for this virtual cell. This merging procedure 
ensures that the area of the control volume is always greater that \AxAy. 

It is convenient to denote by cxf ■ the set of indices {i',]') such that all the virtual cells 
Vti! ji{t) are merged together. For example, if and ^ij+i merge, then = {{i,j), {i,j + 
!)}• 

The previous algorithm dehnes the control volume at time T*. For t > f”, the virtual cells 
change (since the solid boundary moves), and the control volume as well. For t between t'^ 
and the time dependent control volume is dehned as follows: 




u 


rij/j/(f). 


(13) 




In other words, the set of virtual cells selected at time for merging dehnes the control 
volume up to We point out that if the shape of the virtual cells (and hence of the 

control volume) can vary in time, the set cTjj is hxed for t G 

At time we have to take into account that there are new virtual cells, some others 
have disappeared, and the shape of all of them have changed: therefore, a new control 
volume, denoted by has to be constructed (by the previous recursive algorithm). 

We refer to hgure for an illustration of this algorithm. 

While the previous procedure might look complicated, note that most of the virtual cells 
do not merge, and hence C^ j{t) = fiijit) for most of them. 

The area of the control volumes and are computed easily with 




E 

(i'j')eo-r 




and 


77 + 1 ,* 


E 




Note that there is a kind of redundancy with this approach: indeed, in the previous 
example, since and belong to the same control volume, then the control volumes 

and (t) are the same. However, this makes the implementation much simpler, 
while the overhead of the computational time is very small: indeed, the number of merged 
cut cells is very small as compared to the number of gas cells. 
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3.3 Numerical scheme 


From now on, calculations are detailed with the reduced distribution function /. The same 
analysis can be done with g. The cut cell method is based on a finite volume scheme. 
One time iteration (which will be divided into three steps) consists in computing the average 
value /”2 p of the distribution function over the virtual cell from the average value 
dehned by 

= ^ L /p(c,x)ds. 




(14) 




Similarly, and stand for the average values of the distribution function over the 

control volumes and C)(^(t"+^). They are defined by 


= /Pr,x)d^ and /S" = Tm- 


/C" (i-+i) 


U{e,^)ds. (15) 


The hrst step of the method is the computation of p through the average values of / 


over the virtual cells included in Dehnitions (13), (14) and (15) readily lead to 




(16) 




The second step is the time integration of the integral form of the Boltzmann equation ([^ 
between f” and this relation is applied by choosing the surface S{t) as the control 


volume and by using (15). This gives 

n+l,* rn+1,* 


j-n+l 








where the transport and collision terms are dehned by 


(17) 


Tij^p{t) — 


Qi,j,p(t') 


laciAt) 


C?At) 


w(t)) ■ n(t) /p(t,x) dl, 


[Fp(t,x) -/p(t,x))d^. 


(18) 


The transport integral can be computed as follows. First, dehnition (13) implies that the 
integral over dC^j{t) is the sum of the integrals over the contours of all the virtual cells Qiij/{t) 
that merge into the control volume Cf'j{t). Moreover, the velocity w ■ n is zero for the four 
edges that ht with the Cartesian mesh lines, while this velocity w is for the last edge of 
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the cell, since it fits with the solid boundary. Finally, the transport integral is written as: 








(vp - w(t)) ■ n(t)/p(t,x) dl 


(vp - u^(t)) ■ n(t) fp{t, x) d/ + / Vp ■ n^iit) fp{t, x) dl 
(i) JL{t) 


(19) 


where L = ULpj,_i ULp j,_i is the union of the four edges that £t with the 

Cartesian mesh lines. The collision integral of (18) is readily approximated by ~ 


The approximation of the time integralTn the right-hand side of (17) will be detailed 
in the following sections. 


-^n+l 


The third step of the method is the computation of fij^p, the average value of / in the 


new virtual cell ■ This is done by distributing the value given by (17) to the cells 


merged into the control volume 

_pn+l,* 

7 i,j,p ■ ~ J i,j,p ■ 

A first summary of the cut cell method is given below: 






( 20 ) 


1. The virtual cells flijit) merge into some control volumes Cij(t) and the values are 
computed with (16). 


2. The numerical scheme (17) is applied in order to computed the values f^^^p’*- 


“ 7771+1 


3. The values f^ jp are updated with formula (|20|). 


The three steps of the method are illustrated in hgure for various situations. Note that 
because of the merging procedure, there is no issue of appearing/disappearing gas cells: in 
other words, a small virtual cell necessarily merges with a larger cell before it disappears, 
and conversely, the average value of / in a new appearing virtual cell is naturally dehned 


through steps 2 and 3. This ensures that the method is conservative, see section 3.4 


It remains to explain how the time integral is approximated in ( [I^ for step 2, which is 
done in section 3.3.1 and 3.3.2 and to explain how the motion of the solid body is taken into 


account: this is done in section 3.3.3 The complete scheme is summarized in section 3.3.4 


3.3.1 First order explicit scheme 

A backward Euler method is applied in order to get a hrst order approximation of the 


Boltzmann equation (17). This means that the time integral in (17) is approximated by the 
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rectangle rule. Using (19) and (18), we find that relation (17) becomes: 


/, 


n+1,* 


p _ At 


n+1,* -I i,3,p 


E 


n+1,* 


-nn _ -rn 

X! A! m X,! Af ^ 

* T 2 ^ 2’’^ 'P 


~rn 

(f} . A + 

_L 'I 

^ n+1,* n 


V 1 + -Ei' 


( 21 ) 


where *U-^i and are the numerical fluxes across the hve edges of the virtual 

cell that are computed with a standard upwind scheme: 




:= \L 


i+ij-p • I (min(np,, 0) ^ + max(r;p,, 0) , 

+i,p := I4"i+i I [min(r;p„ 0) + max(np„ 0) //^-J , 

■= l^^il [min([vp - u^(r, r"^.)] ■ n^(r, r^^.), 0) /^(r, r”^-, Vp) 
+ max([vp - u^(r, r”j)] ■ n^(r, 0) /^J , 


( 22 ) 


where is the center of L^j. It is recalled that Vp-^ and Vp^ are the coordinates of the 
microscopic velocity, i.e. Vp = {vpj^,Vp^). Moreover, n^(t"',r”^) is the outward normal to the 
edge L” ■, that £t with the physical boundary. The computation of the velocity r” •) of 


the boundary is detailed in section 3.3.3 Finally the discrete boundary condition is similar 


to its continuous form (10), that is to say: 


/^(r,rT e r,Vp e M.) = 0M[l,u^(r,rT),T^], 


RT 


g^{r,rl^ e r, Vp e Un) = (/.^M[l,u^(r,rT),T^], 


(23) 


where 0 is given by 


EvpeVout(’^p “ U”i)/”i,pAv,,Anp 


E 


v^eVinAp 




Note that the boundary condition is defined only for the velocities Vp G Vm = {v|(v — 
u^) ■ n^; < 0}, which is compatible with the dehnition of the numerical boundary flux ^Zi,p 


(see (22)) 


3.3.2 First order semi-implicit scheme 


When the Knudsen number is very small, the previous scheme (21), which is explicit, is too 


expensive: the CFL condition induces a time step which is of the order of the relaxation 
time. It is now standard in kinetic theory to use instead an implicit/explicit scheme (see 
for instance [32] for the BGK equation and [T3| for other methods). The idea is to use an 
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implicit scheme for the stiff collision part, while the transport part is still approximated by 
an explicit scheme. This kind of scheme can be easily extended to the cut cell method. For 
instance, the simplest hrst order semi-implicit scheme is 
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where the numerical fluxes are still computed with (22). The equilibrium function 
depends on the macroscopic quantities that have to be computed before this can be 


done by summing (24) over the discrete velocities, so that the collision term vanishes, which 
gives an explicit relation for these macroscopic quantities. 


3.3.3 Motion of the solid body 


The motion of the solid body is taken into account in the scheme by the variation of the 


area of the control volume (from to and by the velocity of the solid 


boundary, see (21) and 


These quantities are computed as follows. 


Let c{t) and 6{t) be the coordinates of the center of mass and the inclination of the solid 
body. Its translational and rotational velocities are then denoted by c{t) and 0{t). The 
motion of the solid body, with mass m and moment of inertia J, is modeled by the Newton’s 
laws of motion that are discretized as follows: 


cn+i = c" + At c" and 0”+^ = + At r, (25) 

^n+i = and 0^+^ = ^ + At T^/J, (26) 

where F and T are the force and torque exerted by the gas on the solid body. They can be 
computed by using the stress tensor at the boundary with the formula 


F 


E./.ri./; dl 


'dn„ 


and T 



c) X (E^n.^;) dl. 


These relations can be approximated by any quadrature formula, and we hnd it convenient 
to use a summation over all the cells of the computational domain to avoid too many tests. 
This yields: 
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(27) 

(28) 
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Since |L"j| is non-zero only for solid edges of cut cells, these formula are consistent approx¬ 
imations of the previous dehnition. 

Moreover, while the stress tensor is dehned by (12), the boundary condition has to be 
taken into account to dehne the distribution of incoming velocities, and we set 

'■“j) = X) ("'P “ r^)) ® (Vp - rj,, Vp) Av^Avy 


E ( 

VpGVout 


u»(C, r”,)) 0 (Vp - u„((”, fy ))/”,,p Av,Avy. 


(29) 


The boundary condition r”^-, Vp) is dehned by (23). The new velocity of the wall is 

hnally computed with 


u^(r+\ rl+^) = + (r”+^ - (30) 

where a-*- is the vector obtained after a rotation of 90 degrees of any vector a in the counter¬ 
clockwise sense. 


3.3.4 Summary of the numerical scheme 

For the convenience of the reader, the different steps of the complete numerical scheme are 
summarized below. 

We assume that, at time t^, all the following quantities are known: the average value 
of the distribution function ^ ^ in each virtual cell the parameters of position (c”',6*"') 
and velocity (c”, 9 "') of the solid body, and hence the wall velocity u^„(t"', r”^). One time 
iteration of the numerical scheme is decomposed into the following 7 steps: 


1. The position of the solid body at is computed with (25). 


2. The virtual cells are arranged into control volumes C^j{tn) following the rule given in 
The distribution function is averaged over the control volumes 


section 

3.2.3 

see ( 

15 

. 


3. The virtual cells and control volumes are moved according to the new position com¬ 
puted at step 1. The areas sfj and are computed by using a level set method 

(see appendix [A| for more details). 


4. The value of the distribution at each solid boundaries is computed through boundary 


condition (23). 


5. The stress tensor Ti^ at each solid boundaries is computed with (29), which gives force 
F” and torque T" with (27) and (28). Then the translational and rotational velocities 


are computed at time by the discrete Newton laws (26), and hnally the new wall 


velocity is computed with (30). 
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6. Scheme (21) is used to pass from ^ (the average value of / at time in the control 
volume to (the average value of / at time in the control volume 




-^n+1 


7. The values f 'ijp of / at time in each virtual cell that are merged into the control 
volume are updated with (20). 


3.4 Properties of the scheme 

3.4.1 Positivity 


Standard arguments show that the explicit scheme (21) preserves the positivity of the solution 
if At satishes the following CFL condition: 


At I max 
id 
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id,P 
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(32) 


+ \Li',j'\{{^p - u^(r,rT)) ■ n^(r,rT))+) 


and = max(n,0), v~ = min(n,0). For a correct description of the flow, it is necessary 
that the discrete velocity grid contains the solid velocity Ui„ at any time (since the diffuse 
boundary condition produces particles with as mean velocity). Under this assumption, 
it can be proved that (j)ij,p < C{vx,ra!al + Vy^max!^y)) where C = 1/2. This gives a 
simpler CFL condition, and in practice, it is relaxed by taking C = 0.9 without any stability 
problems 

For the semi-implicit scheme (24), a similar condition can be found, which is independent 
of r, 




3.4.2 Conservation 


Since our scheme is a hnite volume method in which a conservative reflexion boundary 
condition is applied to compute the numerical fluxes at each solid edge, it is naturally 


conservative. This is proved below for the explicit scheme (21). 

Let M" be the total mass of gas in the gas domain Vtg 
write the total mass at time as 


at time t'^. It is convenient to 


N^,Ny 

Mn+1 ^ 

id=^ 


N^-1 
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where = 1 if Xjj G fig and 0 else: this function allows to take into account merged cells 
of a same control volume only once. Indeed, there is only one couple of indices in a, 

for which Xj/y is inside the gas domain. 


n 


Then is replaced by its value given by (21) and we get 


N^,Ny 7V2-1 

Mn+1 

*J=1 P=0 

Indeed, opposite fluxes across same Cartesian edges cancel out, the velocity sum of the 
collision operator is zero, and there remains only the numerical fluxes ^ across the solid 
edge of cut cells. By using the boundary condition, the velocity sum of such fluxes gives 

V2-l 

p=0 VpSVin 


E ( 

VpGVout 


rT)) ■ n^(r, rT) /”-p|TT| 


= 0 . 


This shows that = M'^ and concludes the proof. 

Finally, note that it is standard from the positivity and conservation properties to con¬ 
clude that the scheme if stable. 


4 Numerical results 


4.1 Translational motion under radiometric effect 

An inhnite set of thin plates of height D is located in an inhnite channel of width AD. The 
distance between two plate centers is AD. This experiment has been numerically investigated 
in 021 , where the plates are inhnitely thin (their thickness is zero). The temperature of the 
right side of a plate is twice as the temperature of its left side. This induces a force that can 
be interpreted as the difference between the radiometric force and gas friction, and make 
the plates move. We point out that all the plates move at the same velocity. After a while, 
the radiometric force is balanced with gas friction and the total force decreases to zero. At 
this time, the plates have reached their stationary velocity. In |12] the ES-BGK model is 
used and simulations are made in a moving reference frame in which the channel velocity 
is positive while the plates are motionless. The total force applied on the plates is zero for 
a specihc channel velocity which corresponds to the stationary velocity of the plates. In 
our work, the BGK model is used to described the behavior of the gas and plates are D/10 
thick (see Fig. 5(a)). The transient motion of the plates (i.e. plates velocity before the 
stationary state is reached) can be simulated with the cut cell method in a hxed frame of 
reference. In that case, the motion depends on the mass of a plate. If we denote by po fhe 


16 




initial density of the gas, this mass is set to m = po-D^/2. The computational domain is 
a rectangle of size AD x 2D that describes the upper part of the channel (Fig. |5(a)[ ). The 
left and right boundary conditions are periodic so as to simulate the inhnite channel. The 
bottom of the computational domain hts with the center of the channel: specular-reflection 
boundary condition is applied to take into account the symmetry. At last, the top of the 
domain corresponds to the wall of the channel, the standard diffuse boundary condition is 
used. During a whole simulation, there is always one and only one plate in the computational 
domain: when a plate goes out, an other one comes in. Note that solid and cut cells only 
appear near this plate. To take its temperature into account, the boundary of the plate 
is modeled with the diffuse-reflection condition. Finally, the relaxation time and Knudsen 
number are dehned by the relations: 


h 

pRTq 


and Kn 


2 ^/WTq 
y/n poRTo/p 



where p is the viscosity of the gas and Tq the initial temperature of the gas. 

Our simulations have been done for a wide range of Knudsen numbers from 10“^ to 1. 
Converged results are obtained for a velocity grid that contains from 20^ to 40^ points (this 
depends on the Knudsen number) and for a spatial mesh made up of 400 x 200 cells, which 
means that a plate encloses 10 cells. For coarser grids, the number of cells enclosed in the 
plate is too small to capture the shape of its edges with enough accuracy. As expected, 
the magnitude of the velocity of the plates increases until they reach their hnal velocity, as 
illustrated for three different Knudsen number on hgure |5(b) 

The variation of the stationary velocity of the plates is plotted as function of Kn on 
hgure For small Knudsen numbers, the hnal velocity seems to be proportional to a/Ku 
while this velocity tends to a constant for high Knudsen number. These simulations show a 
good agreement between our results and the results obtained in [42] . 


4.2 The Crookes radiometer 

The Crookes radiometer was invented by Crookes in 1874 [12]: it is a glass globe containing 
four vanes immersed in a low pressure gas. Each vane has one black side and one shiny side, 
and when the globe is exposed to light, the vanes rotate. This was hrst understood as a 
rarehed gas dynamics ehect by Reynolds |3l], but there are still discussions on the order of 
magnitudes of the forces involved in this device. We refer to the recent review of Ketsdever 
et ah [21| for historical details. Recently, numerical simulations improved the understanding 
of the radiometric effect, like in mi 1121 EZ]. 

The dynamical acceleration process of the vanes has been recently studied in [9] : it uses 
the unihed gas-kinetic scheme combined with a moving mesh approach [101 to simulate the 
motion of the vanes. In this case, the moving mesh approach is very convenient because the 
initial mesh just rotates without distortions. In this section, we show that the same results 
as [2] can be obtained with our cut cell method. 
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Since this test is for illustrating our 2D method, the device simulated here is a 2D 
radiometer composed of endless vanes immersed in an unbounded cylinder of radius R = 
20cm. The length of a vane is L = 0.1cm, its thickness is I = 0.01cm (figure 7(a)) and 
its rotational moment of inertia is J = 4.9 x 10“®kg-m^. The temperature of the black 
side of the vane is supposed to be higher than the temperature Tc of its shiny side. These 
temperatures are set to Th = 400K and Tc = 350K while the temperature Tq of the globe 
is equal to Tq = 300K. Note that all the boundary conditions are computed with diffuse 
reflection conditions. In order to compare our results to [U], we take the same Shakhov 
relaxation model with Pr = 2/3. Moreover, the relaxation time is computed by the equation 


P \Tr 


(34) 


where /i is determined by the hard sphere model for argon which yields /i = 1.678 x 
10“^Nsm“^ and u = 0.68. Finally, we set the initial density po = 8.582 ■ 10“®kg-m“^ and get 
a Knudsen number based on the length of a vane equal to 0.1. 

Converged results are obtained with a 30^ points velocity grid and 400^ cells spatial mesh. 
This large number of spatial cells is required to describe the shape of the vanes with enough 


accuracy. We plot in hgure 7(b) the radial velocity of the vanes as a function of time and 
compare our results to those obtained in [U]. The results are in very good agreement. 


4.3 Roots blower 

Various kinds of vacuum pumps are used in industrial processes. A common one is the Roots 
blower. It is made of several lobes that rotate simultaneously. As a result, the gas is trapped 
by the lobes at one side and then carried to the other side of the pump. The simplest shape 
of Roots blower is a two-lobed rotor. In this case the prohle of a lobe is dehned by sections 


of epicycloid and hypocycloid (figure 8(a)). In parametric coordinates, this profile is given 
for all 9 G [—vr, — f] U [0, |] by the epicycloidal equation 


x{9) = 5r cos{9) — rsin(56*), 
y{9) = 5rsin(6*) — rsin(56'), 

and for all 6* G [—|, 0] U [|, vr] by the hypocycloidal equation 

x{9) = 3r cos{9) + r sin(36'), 
y{6) = 3rsin(6*) — rsin(36'). 


where r is the radius of the small generating circle that rolls on the large circle of radius 
4r. In our simulation, we took r = 3.8cm for both lobes. The full geometry of the pump is 


detailed in hgure 8(b) Note that lobes are not in contact : the minimum distance between 
them is d = 1.6cm. 
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Since this type of pump mostly operates in atmospheric environment, the initial condi¬ 
tions are given by Tq = 300K, Pq = lO^Pa, Uq = 0, and the considered gas is argon. The 


relaxation time is computed with formula (34) where the viscosity coefficient and index for 
argon are provided by Bird |H], that is /i = 2.117 x 10“^^Nsm“^ and u: = 0.81. Because there 
is no friction between the lobes, a Roots blower can proceed at a rotary speed that range 
from 1500rpm to 3000rpm (~ 150rad-s“^ to 300rad-s“^). For the following simulations, the 
rotational velocity of the lobes is set to 6^ = ±200rad-s“^. Note that since the velocity is 
imposed here, step 5 of the algorithm is not used (see section 3.3.4). In practice, the tem¬ 


perature of the lobes tends to increase because of the mechanical heating due to their high 
rotary speed. However, to make it simpler, the wall of the Roots blowers and its lobes are 
modeled with diffuse boundary conditions with constant temperature Tq. At the left side of 
the pump, it is assumed that all the gas surrounding the computational domain is in the 
same state as the gas located at the inlet. This can be modeled by a Neumann boundary 
condition. At the outlet (right side of the pump), we assume that the gas is released in 
the atmosphere, and hence the boundary condition is given by a Maxwellian built with the 
initial conditions Pq, Tq and uq. 

Pressure contours at several times are shown hgure We observe that the pressure at 


outlet does not change while the pressure at inlet decreases, as it is shown in hgure 9(d) 
We stop the computation at t = 0.1s. At this time, the inlet pressure is 80% of the initial 
pressure, which means that we get a pressure drop of 20%. 

We point out that this simulation is only a qualitative analysis. Kinetic equations are 
not really relevant here because the Knudsen number is very small: Kn 3 x 10“® for 
a reference length equal to the distance d between the two lobes (it would be larger with 
a smaller distance). Hence Navier-Stokes equations might be more relevant in this case. 
However, this simulation shows that the cut cell method works well with complex shaped 
objects for moderate velocities hows (Ma ~ 0.15), while this would be more more difficult 
with the moving mesh approach, for instance. 


5 Three dimensional flow simulations 

In this section, the cut cell method presented in section is extended to 3D simulations. The 
approach is only detailed here for the simulation of the Crookes radiometer, which means 
that only pure rotation is considered. The angle between the position of a vane at time t 
and its initial position is denoted by 6{t), and 6{t) stand for its rotational velocity. The gas 


governing equations are that of section 2T and the velocity discretization is similar to the 
one explain in section 3T there are velocity points and the velocity is denoted by 
Vp = (up^, Wps)- The different steps of the cut cell method - that consis ts in updating the 
values 6*"', 6*” and p for the next time step - are that of section 
below to highlight the differences with respect to the 2D case. 

We assume that at time we have the average value of the distribution function F^ -^^ p 
in every virtual cell Like in 2D, a virtual cell is dehned as the intersection between 


3.3.4 


They are written 
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the gaseous domain and the cuboid cell of the Cartesian mesh, but now it is a 

polyhedral with seven virtual faces. The hrst six faces £t with the Cartesian mesh interfaces 
and are denoted by S'” i .,, S'”., i ,, , , i. The seventh one is a plane approximation of 

2J±2,A:’ *J,A:±2 ^ 

the solid boundary denoted by S'jj^fc and is its normal vector directed outward. We also 
assume that at time f”, the angle 0” and rotational velocity 0” of the solid are known. The 
different steps are the followings: 


1 . 

2 . 


The position of the solid body at is computed with (25). 


The virtual cells are arranged into control volumes 
in section 


following the rule given 


3.2.3 


which is naturally extended to 3D. The average value F[^j^k,p f^e 
distribution over the control volume is computed: 2D relation (16) is replaced 

by 
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F 






E 




T/ r 

^ J',k' ^ i',j' 


i,j,k 


3. The virtual cells and control volumes are moved according to the new position com- 

puted at step 1. The volumes of the virtual cells control 

volumes and at times f” and respectively, are computed. 

4. The boundary condition is computed for each cut cell with the discrete form of ([^, 
which yields: 

FwiF.Fljk e r,Up e Vm) = 0A4[1 ,m^,T^], 

where (j) is given by 


• nwiF:^,j,k)Fl^j,k,p^Vx^Vy^Vz 

Here, ff^k is the center of Sij^k and Vin = {vp\{vp - a,)) • < 0} 

and Vout = {vp\{vp — u^iF^Fljk)) ' ^w{FiFljk) > 0} are the sets of incoming and 
outgoing velocities, respectively. We introduce the cylindrical coordinates of r^j k = 
(r” cos a”, r” sin a”, z^) in order to write the boundary velocity as Uw{F, ^ 

(sin a”, cos a”, 0). 


5. The stress tensor is computed by 3D extension of (29): 

fc) = (fV - f^j^k)) ® (vp - Flj k))F^{r, fl - k, Vp) Av^AvyAv^ 

Vp&Vin 

+ ^ ^ {Vp — Uyj{t ,PiJ^k)) ^ {"^P ~ ,r^J^k))FiJ^k,p ^'^xAVyAVz- 

l^pGVout 

The rotational velocity is then computed with 

. , N, _ 

j,.+i =r+—Y,Y.T. Hi* ^ 


i=l j=l k=l 



1 

o 

_1 


0 

1 


20 











where the sum is an approximation of the torque acting on the vanes. Then the new 
wall velocity at time is computed. 


6. Scheme (21) is easily extended to 3D to compute the average value of / at time 

fn+i control volume indeed, we write the integral form of Boltzmann 

equation (j^ with V{t) = and simplify the transport term just like in (19) to 

get the hrst order explicit scheme 
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(35) 


where the upwind numerical flux i .. , . i , , . i and 

1-\- 2 p i^iP p iP 


ij,k,p 


are: 
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[min(np^, 0) + max(np^, 0) 

[min(np,, 0) Tl”+i,fc,p + max(np^, 0) J 

[min(np^, 0) Tl^-fc+i^p + max(np^, 0) 

[min([np - 0) ^ Vp 

[([up - fc)] ■ 0) T;". ^ J 


max 


7 =-n+l 


7. This average value of the distribution function is distributed to the virtual cells 
merged into Cl^p{t'^+^) by F'^^Kv ■= Ktk^' 


5.1 Octree procedure 

Full 3D simulations are computationally very expensive. For instance, based on the 2D 


computations presented in section 4.2, a 3D simulation of the Crookes radiometer requires 
200^ degrees of freedom for the space discretization and 30^ for the velocity discretization. In 
addition to these 2 160 x 10® grid points, 10 000 time steps are expected to reach the stationary 
rotational velocity of the vanes. In conclusion, even a massively parallel computing is not 
sufficient to do to this simulation in a reasonable computational time (i.e. less than several 
weeks). 

In order to reduce the computational time, an octree procedure is implemented. First, a 
coarse Cartesian mesh is initialized and rehned around boundaries. This means that all the 
cells that are close enough to the boundary are divided in 8 smaller cells (or 4 in 2D). These 
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new cells can be divided again if necessary; and we define the depth of a cell as the number 
of divisions that lead to this cell. The splitting criterion is the following: 


if 


<2V2^ 

2 ^ 


and d < dmax the cell is divided, 


(36) 


where 0 is the distance from the cell center to the nearest boundary and dmax is a prescribed 
maximum depth. For the 3D simulation of the Crookes radiometer, the coarse mesh is made 
of 50^ cells and d m ax = 2. With these parameters, the spatial mesh only contains 240 000 
cells. Since the boundary moves, the mesh is adapted to the new location of the boundary 
at the beginning of each time iteration. Note that it may happen that 8 cells have to merge 


during this process, if all of them no longer respect criterion (36). 

When the mesh changes - i.e when cells split or merge - the distribution function has to 
be interpolated on the new mesh. This is done by assuming that the distribution function is 
constant over a cell and by using a standard restriction/prolongation method (by average and 


0th order interpolation). Since every cells are cuboid, scheme (35) can be applied by using 
an appropriate data structure to access the neighboring cells of each numerical interfaces: 
here, we use the standard Z-ordering which is very efficient for that. 


5.2 Parallel implementation 

We describe here two natural strategies for a parallel implementation of a kinetic solver with 
the Message Passing Interface (MPI) library, and we propose our own hybrid technique. 

The first method is velocity parallelization, or decomposition domain method in the 
velocity space: each processor computes the distribution function in the whole space domain, 
but for only a part of the discrete velocities. This approach is for instance used in [1] and 
more recently in |l6]. For a given discrete velocity, the scheme is independent of the other 
velocities: then it is is used independently by each processor, and each of them compute 
partial moments (by using its own reduced set of discrete velocities). To compute the full 
moments, all processors gather the sum of the partial moments. 

The second method is a more standard space domain decomposition, see [SUES]. Each 
processor computes the distribution function for the whole discrete velocity grid, but only 
for a subdomain in the position space. To compute the numerical fluxes across the interfaces 
between different subdomains, the method requires communications: it is sufficient that each 
processor sends the distribution of its interface cells to its neighbors. We refer to |11] for a 
comparison of these two strategies. 

Our implementation combines these ideas. Each processor uses the scheme on a space 
subdomain, for a partial set of discrete velocities. The space domain decomposition is made 
on the initial mesh, before the refinement procedure: the Cartesian structure of this mesh 
makes the portioning very easy. For problems with moving boundaries, it is difficult to 
ensure a good dynamic load balancing between different processors: the number of cells of 
a subdomain can change a lot due to the space refinement induced by the displacement of 
the solid obstacle. In order to optimize the workload distribution, we use a small number of 
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subdomains. Groups of processors are given to each subdomain, and each processor will apply 
the scheme for a partial set of discrete velocities. The advantage of this technique is that 
we can use a large number of processors without a two large number of space subdomains. 
This makes the workload well balanced during the simulation for each subdomain, at least 
for the Crookes radiometer presented in the following section, since there is always one vane 
in each subdomain. For the corresponding 3D simulation, our technique is quite efficient, 
since it has been made with 240 processors for a CPU time lower than 12 hours. 

Note that there is an other kind of hybrid parallelization which uses both MPI and 
OpenMP libraries (see i), but this is not what is used here. 


5.3 Numerical example : the Crookes radiometer 

First, the implementation of the method has been checked with a 3D simulation of the plane 
2D radiometer similar to the one presented in section 4.2[ this 2D geometry is extruded 
to get a cylinder shaped radiometer (see hgure 10, left), and periodic boundary conditions 


are imposed at the upper and lower boundaries to simulate the inhnite vanes. In this case, 
the moment of inertia is 97po-h^/370, where L is the height of the extruded vanes. For a 
Knudsen number of 0.5, 2D and 3D simulations give exactly the same results, as it is shown 
in hgure right. For this comparison, the AMR technique is used for both simulations. 


and a plane section of the 3D mesh is the same as the 2D mesh. 

Now the real 3D radiometer is made of four square shaped vanes. The dimensions of 
the vanes are L for the diagonal and L/10 for the width. They are immersed in a sphere of 
radius 2L, and their centers are in the plane z = 0, at a distance 0.75L from the center of 
the sphere. The corresponding geometry is shown in hgure The moment of inertia of the 
vanes is computed with a material density po equal to the mass density of the surrounding 
gas which gives J = -^ PoL^- This density is not realistic, but it makes the vanes faster, 
and it is easier to observe their movement. 

At t = 0, the radial velocity of the vanes is zero and the temperature Tq in the domain 
is uniform. The gas is governed by the BGK model, where the relaxation time is r = 
boundary conditions are diffuse rehections with constant temperatures: 
the sphere boundary is maintained at temperature Tq and the white side of the vanes as 
well, while their black side is maintained at temperature 2To. On the edge of the vanes, the 
temperature is discontinuous (Tq on one part and 2To on the other part). 

We also use the 2D simulation to estimate the resolution required by the 3D computation. 
The difference between the results obtained with a 2D mesh of 50^ cells rehned by the AMR 
technique with a maximum depth dmax = 2 and the hne Cartesian structured mesh of 500^ 
cells is less than 5%, which is considered as sufficiently accurate here. Consequently, a 3D 
AMR mesh of 50^ cells with a depth of 2 should be accurate enough for a 3D simulation. This 
is computationally possible, since this mesh contains 240 000 cells, which is much smaller 
than the equivalent Cartesian mesh of 200^ = 8 000 000 cells. The corresponding simulation 
is shown in hgure at different times. 

We have made other simulations with three Knudsen numbers 0.1, 0.3, and 0.5 until 
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the steady state is reached. In hgure 13(a) is shown the evolution of the radial velocity: we 
clearly see the convergence to a constant velocity, which is larger for Kn = 0.5. A comparison 
between 2D an 3D geometries is shown in hgure |13(b) the 3D vanes are clearly faster that 
the 2D vanes. 

Note that this test is just shown to illustrate the potential of our method. We are not 
aware of any similar simulation in the literature so far, and we are not able to present any 
comparison. Moreover, we do not claim this is a realistic simulation, since the moment of 
inertia of the vanes is too small, and their width is too large. However, we do not know 
any experimental measures of the motion of the radiometer. If any, it would probably be 
necessary to make a more intensive simulation, since the rehnement should be stronger 
around the vanes that are generally very thin. 


6 Conclusion 

A numerical method for solving kinetic equations with moving obstacles has been presented. 
This method is an extension to the kinetic theory of the cut cell technique used in com¬ 
putational fluid dynamics. The main advantage of this algorithm is that it combines the 
simplicity of the Cartesian grid based methods to the accuracy of the body htted methods, 
which ensures exact mass conservation. The method is easily extended to 3D flows, and its 
accuracy has been proved with the simulation of a Crookes radiometer. Another advantage 
of our approach is a simple and generic treatment of all kinds of cut cells. This is essential, 
especially for 3D problems in which there are many different kinds of cut cells. 

Our goal is now to improve the accuracy of our method by using a second order scheme. 
Since the mesh is Cartesian, the main difficulty is to approximate the gradient of the distri¬ 
bution function on the cut cells with enough accuracy. Such an extension has already been 
done (see [I3]), but while it works well for non moving obstacles, it is not efficient enough for 
general problems. An other perspective is the validation of the 3D algorithm, in particular 
by using experimental data. 
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A Computation of the cell geometric parameters 


We describe below how these parameters are computed for 2D problems. There is no specihc 
difficulty to extend our algorithms to 3D problems: however, the formula are a bit long, and 
to shorten the paper, this extension is left to the reader. 

The computation is made with the following four steps: 

1. identihcation of the type of each cell (solid, gas, our cut cell); 

2. for each virtual cell, computation of the lenghts of its hve edges, and computation of 
its normal vector of its hfth edge (the one which is a linear approximation of the solid 
boundary); 


3. identihcation of the cut cells that have to merge, and computation of the set aij] 

4. computation of the area of each virtual cells and corresponding control volumes. 


Before we describe these steps in the following sections, we point out that a level-set 
signed distance function 0 : x —>■ 0(x) is systematically used. It gives the shortest distance 
between a point x in the computational domain and the solid boundary. This function is 
negative if x is inside the solid, and positive if it is inside the gas, and hence the zero level-set 
of 0 is the solid boundary. In our algorithm, the values of 0 are computed analytically at 




of the Cartesian grid, and we set := 0(xjj_ij_|_i). Since solid 


each node 

boundaries move, it is necessary to update these values at each time step. 

From now on, we consider a single cell Djj, and in order to simplify the notations, indices 
i and j of variables dehned at the vertices of this cell will be omitted. For instance, the values 
and will be denoted by x+^+ and 0+,+ . All the notations used here are 


^i+k,j+h 


shown in hgure 14 


A.l Identification of cell types 

It is clear that flij is a gas cell if the values of 0 at its four vertices are positive. At the 
contrary, it is a solid cell if these values are all negative. Finally, if these values have different 
signs, the cell is cut by the solid boundary. These three types are identihed by looking at 
the sign of m := min(0_^_, 0+_, 0+,+, 4>-,+) and M := max(0_ _, (j)+-, 0+,+ , 4>-,+)- 

m > 0 -v^ flij is a gas cell, 

M < 0 -v^ Djj is a solid cell, 

M > 0 and m < 0 Djj is a cut cell. 

Of course, if a cell is crossed by a solid object thinner than Ax, it will be considered as a gas 
cell: it is important that the Cartesian grid is hne enough to resolve all the solid objects. 

Finally, the case = 0 is complex and is avoided by using the modihcation 0 : = 
sign(0) X max(| 0 |,10“^°Ax). This means that grid points that are exactly on the solid 
boundary are numerically considered as moved on a distance of Ax. This has no 
influence on the accuracy of the results. 
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A.2 Lengths of edges and normal vector of the virtual cell 

To each of the four edges of are associated the four points et dehned as follows. 
If it is a cut cell, two of these points are intersection points of an edge with the solid boundary, 
and the two others are not used by the algorithm. When hljj is a gas or solid cell, none of 
these four points is used. These points are dehned by linear approximations 





^.4 

1 

> 

1 

H- 

y "' — 

V 1 Av 

and x))^ = 

0+,± - 0-,± 


(/o-i 1 1 

" 0±,+ - 0±-. 





To compute the lengths of the edges of hljj, we use the symbol that is 1 if x± ± is 
inside the gas, and 0 if it is inside the solid. This value is given by = max((;/)± O)/|0±^±|. 

Also note that looking at the products (or tells us if left and right edges 

(respectively, upper and lower edges) are crossed by the solid boundary. This notation allows 
us to easily write the four Cartesian edge lengths of the virtual cell hljj: 


IT" 1 .1 

|T"i| 


<^±,-11x1^ - x±,-ll + <^±,+ l|x±,+ - x^jll + (5±,_(5±,+ i|x±,+ - x±,_i|, 
- x_,±|| + (5+,±||x+,± - x^ill + (5_,±(5+,±||x+,± - x_,±||. 


The length | of the hfth edge (that ht with the solid boundary) is the norm of the vector 
dehned by the two solid boundary/edge intersection points. 

Finally, the normal vector to this edge is computed by a Green formula, which naturally 
gives the correct outward direction: 



A.3 Computation of aij 

We consider one virtual cell hljj. The set Uj ^ collects the indices of all the virtual cells that 
merge with Qij. In this set, let ns denote by the indices of the unique virtual cell 

whose center Xj/j/ is inside the gas, that is to say 

^( 0 -,- + 0 +,- + 0 +,+ + 0 -,+) > 0 - (37) 
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We point out that the rule given in section 3.2.3| ensures that this ’’master” virtual cell is 
unique. Its indices {i'ij') are determined by the following algorithm 


while ( |3^ is false, do 


M 

:= max(|L”^i^^.|, |. 


m-^AAL"- 

' *J+2' ' 

if 

\r^ 1 

= M then 


:= (^' + 1,/) 

if 

m 1 -1 

*-2j' 

= M then 

A',f) 

:= 

if 


= M then 


:= + 

if 

\L^. i| 
*J-2' 

= M then 


:= 


end while 


(38) 


Finally, note that the set aij is not really computed: practically, we compute the numer¬ 
ical fluxes across the edges of each virtual cells, and these fluxes are directly added to the 
fluxes of the master cell of indices [i',j'). 


A.4 Virtual cell and control volume areas 

The area of a virtual cell can be computed with the length of its edges and the coordinates 
of its vertices by a Green formula: 




d5= - 




2 /n. 


V ■ xdS = (i+,+x+,+ + 5 +,-X+,_ + yjx”,) . 

+ \\LUJ . 

+ I (^+,+x+,+ + 5-,+x_,+ -1- (5+;+x^+) ■ 

+ - + '^+-x+- + <5^;_x^_) ■ 


0 

1 

0 

0 

1 

0 

1 


+ -Am mmm +^ 


Here, we used 5^ + := 1 — 6±^+6±- which is 1 if left or right edges are crossed by the solid 
boundary, and 0 else. This formula is nothing but the sum on each edge of the dot product 
between its normal vector and the vector pointing to the center of the edge, multiplied by 
the length edge. 
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Finally, the area of the control volume is computed and stored in the master cell of indices 
(F, j') dehned by algorithm (38), with the following loop along all the cut cells: 


Sij = 0 for every cells. 
For all {i,j) do : 


compute with algorithm (38) 




* J 




end do 





Figure 1: Cells of the computational domain are classified in three categories: gas cells 
represented in white, solid cells are shaded and cut cells are hatched. 
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Figure 2: Three examples of virtual cells: the cell 12*^ is drawn with the dashed line, while 

the corresponding virtual cell is drawn with the solid line. This virtual cell is a polygon 

with at most hve edges. Left: hve edges. Middle: three edges, while two virtual edges T” _i 

2 

and T"_i . have zero length. Right: four edges while the virtual edge i . has zero length. 
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Figure 3: Top-Left: the virtual cell (hatched) has to merge with because \ > 

\L‘^_i ■\- Top-Right: this results in the control volume in blue. Middle-Left: at time 

the solid boundary has slightly turned counter clock-wise and the virtual cell is 
smaller. Middle-Right: the control volume at this time now is in blue. Bottom: 

since the larger edge of the virtual cell now is .\^ it has to merge with and 

this results in the new control volume which is different from . 
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Iteration n 


Iteration n + 1 



(a) The boundary is moving from the right to the left. 



(b) Appearing cut cell: the boundary is moving from the right to the left. 
Iteration n ' Iteration n + 1 



(c) Disappearing cut cell: the boundary is moving from the left to the right. 

Figure 4: Illustration of the three steps of the cut cell method. The cells Qij are drawn with 
the dashed line and the virtnal cells with the solid line. The control volnme Cij{t) is 

shaded. 
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(a) Experimental set 

up of the moving 

plates. The computational domain is drawn 
with a doted line. 



(b) Evolution of the velocity of the plates for three 
Knudsen numbers. 


Figure 5: Translational plates under the radiometric effect 
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Figure 6: Translational plates under the radiometric effect: stationary velocity of the plates 
as a function of the Knudsen number, comparison with Taguchi et ah |42j . 
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(a) Experimental set up, the computa¬ 
tional domain is drawn with a doted line. 



Figure 7: The 2D Crookes radiometer. 


35 





















0.4 



(a) Geometry of a lobe. The plain line is 
an epicycloidal and the dashed line is an 
hypocycloidal. 



(b) Pump geometry. All the units are given 
in meter. The computational domain is the 
square [—0.4,0.4]^. 


Figure 8: Roots blower. 
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(c) Pressure distribution at t = 0.09s. 



(d) Inlet pressure as a function of time: Pjniet = 
r ^dx 

J inlet Po 


Figure 9: Pressure in the pump at several times. 
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Figure 10: 3D extruded radiometer: temperature field T/Tq and mesh (left) and radial 
velocity profile as a function of time for 2D and 3D simulations (right). 
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Figure 11: 3D radiometer: cross section in the plane xOy (left) and in the plane xOz (right). 
See hgure for a 3D view. 
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Figure 12: 3D Crookes radiometer for Kn = 0.5: the temperature field T/Tq on the plane 
z = 0 is shown at times t x L/a/ZRTq = 0, 5, 10 and 15 (from left to right and from top to 
bottom). The mesh is shown at t = 0 (top-left). 
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(a) Time evolution of the radial velocities. (b) Stationnary velocity for 2D and 3D simula¬ 

tions for three Knudsen numbers. 

Figure 13: Radial velocity profiles for 3D simulations with 3 different Knudsen numbers. 
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0 - 
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Figure 14: Summary of the notations used in Appendix the cell hlij is shown with a 
dotted line, its corresponding virtual cell flij with a solid line, and the solid boundary with 
a double line. 
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